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We study dynamic features of a plasma consisting of gluons whose infrared dynamics is improved 
by the Gribov-Zwanziger quantization. This approach embodies essential features of color confine¬ 
ment which set the plasma apart from conventional quasiparticle systems in several aspects. Our 
study focusses on a boost-invariant expansion for in- and out-of-equilibrium settings within the re¬ 
laxation time approximation, which at late times can be characterized by the sound velocity, Cs, and 
the shear, p, and bulk, viscosities. We obtain explicit expressions for the transport coefficients p 
and and check that they are consistent with the numerical solutions of the kinetic equation. At 
high temperature, keeping both the Gribov parameter and the relaxation time constant, we find a 
scaling ^/rj oc 1/3 — which manifests strong breaking of conformal symmetry in contrast to the 
case of weakly coupled plasmas. 


I. INTRODUCTION 

The fate of hadronic matter at high temperature is still 
under active investigation nearly five decades after the 
theory of quantum chromodynamics (QCD) was estab¬ 
lished. Above a certain (pseudo-)critical temperature, 
the fundamental quark and gluon constituents are ex¬ 
pected to become the active degrees of freedom, forming 
a quark-gluon plasma (QGP). At asymptotically high en¬ 
ergies, such a plasma has been for long regarded as a 
weakly interacting gas. However, this naive picture has 
been challenged by experimental data [1-6] from ultra- 
relativistic heavy-ion collisions that give evidence for the 
creation of a strongly interacting and correlated system 
in the accessible energy regime [7, 8]. This has raised 
the interest in describing the dynamics of the collisions 
in terms of dissipative fluid dynamics, characterized by 
various transport coefficients [9, 10], see also [11]. These, 
in turn, are sensistive to the long-wavelength features of 
the underlying microscopic theory by force of the Green- 
Kubo relations [12]. 

While the thermodynamic properties of QGD are well 
established at high temperature, both within resummed 
perturbation theory, see [13-15] for reviews and [16-19] 
for the most up-to-date results on thermodynamic quan¬ 
tities, and by lattice methods [20, 21], understanding the 
dynamic properties remains challenging [22] . Let us point 
out the two main reasons. In the small-coupling regime. 
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<7 <C 1, the degrees of freedom successfully described by 
resummed perturbation theory are dressed quasiparticles 
with an effective mass ~ gT [23-26] , where g is the QCD 
coupling constant. The perturbative QCD predictions 
are however in tension with experimental results since the 
expected running-coupling strength in the phenomeno¬ 
logically relevant regime is of the order of unity, g ^ 0{1) 
[13-15]. Secondly, since the initial conditions of ultrarel- 
ativistic heavy-ion collisions are highly anisotropic in mo¬ 
mentum space, the incipient evolution of the created sys¬ 
tem towards a state complying with a hydrodynamic de¬ 
scription is still debated and new formulations of hydro¬ 
dynamics are constructed in this context [27-29]. These 
outstanding issues have in turn motivated the application 
of strong-coupling techniques, based on the gauge-gravity 
duality, see e.g. [30]. One of the striking hallmarks of the 
latter approaches is the lack of distinct quasiparticle ex¬ 
citations in the spectrum of the theory, but rather the 
dominance of hydrodynamic sound modes [31-33]. 

In our present work, we will focus on two main trans¬ 
port coefficients for soft modes in a charge-free plasma, 
namely the shear, rj, and bulk, Q viscosities. In the last 
years, substantial knowledge has been gained regarding 
their properties. In particular, the shear and bulk vis¬ 
cosities were calculated in the high-T regime using per¬ 
turbation theory in [34, 35] and [26, 36], respectively. 
In the phenomenologically relevant temperature regime, 
T ~ (I—4)Tc [14, 15], lattice calculations are still plagued 
by large uncertainties, see [37] for early results on bulk 
and [38] for shear viscosities, respectively. One note¬ 
worthy feature is the enhancement of the bulk viscosity 
close to the critical temperature, which was also found 
in [39, 40]. Encouragingly, the transport coefficients of 
the Yang-Mills theory have been recently computed us¬ 
ing functional renormalization techniques [41, 42]. In a 
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complementary effort, the transport properties of super- 
symmetric plasmas at high temperatures have recently 
been studied using the gauge-gravity duality. For a con¬ 
formal system a lower bound on the shear viscosity to 
entropy density ratio, 77 /s > l/47r, was found [43]. How¬ 
ever, breaking of conformal invariance is necessary for the 
generation of bulk viscosity [44, 45] , and several models 
have been proposed in order to mimic lattice QCD fea¬ 
tures [46, 47]. 

Due to confinement effects, the infrared (IR) regime of 
QCD is strongly coupled, opposite to the case of quantum 
electrodynamics. Long-range correlations in the system 
contribute to the confinement of colored degrees of free¬ 
dom. These features are genuinely non-perturbative and 
are therefore beyond the scope of conventional pertur¬ 
bative techniques [13-15], especially in the study of the 
QGP [48, 49] and for heavy-ion phenomenology. 

Motivated by these considerations, in this work we 
continue our investigation of equilibrium and non¬ 
equilibrium properties of a plasma consisting of gluons 
obtained from the Gribov-Zwanziger (GZ) quantization 
which was introduced in Gribov’s seminal paper [50] 
and later systemized by Zwanziger in [51], for reviews 
see [52, 53]. This prescription improves the IR behav¬ 
ior of Yang-Mills theory by fixing residual gauge trans¬ 
formations that remain after applying the perturbative 
Faddeev-Popov procedure. Consequently, a new scale 7 g, 
the Gribov parameter, is introduced, which leads to an 
IR-improved dispersion relation for gluons [50]. In the 
Goulomb gauge this dispersion relation reads 

E{k) = ^k^ + ^, (I) 

where k is the three-momentum and E is the energy. It 
embodies the expected behavior: a large energy cost as¬ 
sociated with soft gluons, which subsequently amounts 
to the reduction of the physical state space [50, 51], see 
also [54]. The Gribov parameter is an intrinsic Yang- 
Mills scale that is bootstrapped by a self-consistent gap 
equation [50, 51] and explicitly breaks the conformal sym¬ 
metry of the theory. 

This framework has recently attracted a lot of atten¬ 
tion in the theory community. In particular, an excellent 
agreement between Gribov’s result and lattice data in the 
Coulomb gauge has been established for the gluon propa¬ 
gator [55]. As one of the most important features, it was 
early realized that the positivity of the gluon spectral 
function is violated in this theory at any temperature, 
see [56, 57] for reviews. This implies the absence of con¬ 
fined particles in the physical asymptotic spectrum and 
is thus a crucial measure of color confinement. 

The Gribov approach was generalized to finite tem¬ 
perature [58] and it was shown that the Gribov pa¬ 
rameter provides a natural QCD scale that significantly 
improves the IR behavior of the theory and provides 
good agreement with lattice results on thermodynamic 
quantities [59]. Owing to its intimate relation to the 
(chromo)magnetic scale [59, 60], collective massless de¬ 


grees of freedom are present at finite temperature [61], 
see also [62, 63]. Recently, the connection between the 
Gribov parameter and the topological structure of the 
QGD vacuum has been explored, shedding light on a 
profound relationship between the GZ quantization and 
the confinement/deconfinement phase transition [64], see 
also [65]. The impact on real-time observables has also 
been examined lately [ 66 ]. 

Motivated by the results discussed above, the Gribov 
dispersion relation thus provides a unique and straight¬ 
forward way to study the impact of residual confinement 
effects on QGP transport properties. It was for the first 
time employed in the kinetic and hydrodynamic calcu¬ 
lations in our previous work [67], where the analyzed 
system was dubbed the Gribov-Zwanziger (GZ) plasma. 
In addition to obtaining a qualitatively good agreement 
with the Yang-Mills lattice data [ 68 ] in equilibrium, we 
calculated the bulk viscosity in a boost-invariant setup 
[69] using the relaxation-time approximation. The latter 
quantity attracts growing attention in the literature due 
to possible importance of bulk viscous effects in relativis¬ 
tic heavy-ion collisions [70-78]. 

In this paper, we give a detailed discussion of the re¬ 
sults obtained in [67] and, in addition, we present several 
new and complementary results. We study the approach 
to equilibrium of the GZ plasma by varying the charac¬ 
teristic relaxation time, Trei, and demonstrate that the 
late-time behavior of the system can be quantified us¬ 
ing transport coefficients. In particular, we calculate the 
bulk and shear viscosities of the GZ plasma and derive 
their low- and high-T limits. 

One of our important results is the linear scaling rela¬ 
tion C /?7 — 7 tGz(V 3 ~Cs) for the GZ plasma. This relation 
holds for constant Gribov parameter and constant relax¬ 
ation time, and is characterized by the universal values 
Kgz = 5 when T —>■ 0 and Kgz = 5/2 for T —>■ 00 , that are 
independent of 70 - In the intermediate, phenomenolog¬ 
ically relevant temperature range the scaling can be ap¬ 
proximated by a slowly varying coefficient Kgz- A similar 
scaling law was observed from the gauge-gravity duality 
[44, 45] , and differs from the expected high-T scaling ob¬ 
tained within perturbative QGD [36]. 

Our results are compared with the formulas character¬ 
izing a plasma of massive particles obeying Bose-Einstein 
(BE) statistics. It turns out that certain physical observ¬ 
ables (such as, for example, the speed of sound) have an 
apparently similar behavior for the GZ plasma and the 
BE massive plasma, provided the value •\/27 g is adopted 
for the particle mass. However, simple relations do not 
hold in most cases. The observed differences reflect the 
modified IR dynamics that sets the GZ plasma apart 
from conventional quasiparticle models. 

The structure of the paper is as follows. In Sec. H 
we discuss the implementation of Lorentz covariance and 
boost-invariance in our model. The implementation of 
covariance into the model is a fundamental problem con¬ 
nected with the use of the non-covariant dispersion re¬ 
lation, derived by Gribov in the Goulomb gauge, as a 
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starting point. On the other hand, the implementation 
of boost invariance [69] is a simplifying assumption that 
facilitates the analysis of a non-equilibrium system. In 
Sec. Ill, using the concepts of the underlying kinetic the¬ 
ory, we discuss the form of the hydrodynamic equations 
for boost-invariant and transversally homogeneous sys¬ 
tems, denoted below as (0-f 1)D systems. In Sec. IV we 
present the kinetic equation in the relaxation time ap¬ 
proximation and show its predictions for various physical 
observables. The forms of the bulk and shear viscosity 
coefficients are derived in Sec. V, where we discuss also 
the ratio of the two viscosities and relation of C^jr] to the 
speed of sound. In Sec. VD we also comment on the ef¬ 
fect of the momentum-dependent relaxation time on our 
results. We summarize and conclude in Sec. VI. The pa¬ 
per is closed with a series of appendices discussing the 
implementation of Lorentz covariance, various auxiliary 
functions and the low- and high-temperature expansions 
for thermodynamic quantities, the speed of sound, shear 
and bulk viscosities and their ratio. 

Throughout the paper we use the natural system of 
units with c = Ti = ks = The three-vectors are de¬ 
noted by the bold font, the four-vectors are in standard 
font, the dot denotes the scalar product with the metric 
= diag(-|-l, —1, —1, —1). The four-vector defining the 
heat-bath, or the fluid element’s local rest frame (LRF), 
is denoted by u. In this frame = (1, 0, 0,0). 

II. IMPLEMENTATION OF COVARIANCE 
AND BOOST-INVARIANCE 


E{k) appearing in the Gribov-Zwanziger formalism [50]. 
The latter is the energy of the interacting particle with 
a three-momentum k in the reference frame of the heat 
bath. 


B. Thermodynamic-like quantities 

In LRF, the energy density of the GZ plasma described 
by the non-equilibrium momentum distribution function 
f{k) is described by the formula [58] 

^ = 90 j /(fc): ( 3 ) 

where go = 2{N^ — 1) is the degeneracy factor for gluons 
with Nc colors {go = 16 for SU(3)). Generalizing the 
result of [58], the energy density of the fluid is given in 
the covariant form by the expression [67] 

e = y dK E{k ■ u) f{x, k), (4) 

where f{x, k) is the phase space distribution function and 
the integration measure dK is defined as 

JdK{...)^goJ ( 5 ) 

Here, the dispersion relation E{k ■ u) is given by Eq. (2). 
In the case of local thermal equilibrium, the gluon distri¬ 
bution has the Bose-Einstein form [58] 


A. Dispersion relation in covariant form 


1 

exp [E{k ■ u)/T{x)] — 1 ’ 


( 6 ) 


The Gribov dispersion relation, Eq. (1), may be cast 
into the following covariant form [67] 

where u is the four-velocity of the fluid element. We 
assume that the Goulomb gauge as well as the in¬ 
medium value of 7 g are fixed in LRF. We introduce 
= jfcj, which is the magnitude of the three-vector 

k = (fca;, fcy, fc||), and kj_ = yjk"^ + ky such that the re¬ 
sulting four-vector k^ = {k^,k) has standard Lorentz 
transformation properties with k^ = 0. 

It is important to note that the original dispersion re¬ 
lation in Eq. (1) is derived in the Goulomb gauge which 
explicitly breaks Lorentz invariance. In order to regain a 
covariant formalism, which is necessary for description of 
relativistic fluids, one has to make certain assumptions 
regarding the Lorentz transformation properties of the 
quantities appearing in Eq. (1). This problem is pre¬ 
sented and discussed in more detail in Appendix A. Our 
choice, based on Eq. (2), is the most natural and free from 
mathematical ambiguities. We note that in this case the 
zeroth component of the momentum, k^, is different from 


where the temperature T can depend on space and time. 

We proceed in a similiar fashion with the pressure, 
which in LRF is given by [58] 



d^k 


dE{k) 


m 


d?k k^ 
(27r)3 E{k) 



In the covariant version, the pressure reads [67] 



dK 


{k ■ 

E{k ■ u) 



7g 

{k ■ u)4_ 


f{x,k). 


( 7 ) 

( 8 ) 


( 9 ) 


The expressions for the energy density and pressure can 
be used to find the interaction measure (sometimes re¬ 
ferred to as the trace anomaly) 


I = e-3P = 2 


dK 


it 


{k ■ uY E{k ■ u) 


f{x,k). (10) 


This quantity is closely related to the non-conformality 
of the system under consideration. We note that / van¬ 
ishes identically for 7 g 0 (or, more generally, in the 
conformal limit). 
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FIG. 1. (Color online) The temperature dependence of the 
equilibrium energy density e and pressure P (multiplied by 
a factor of 3) for the GZ plasma (red lines) and an ideal 
plasma of massive particles obeying BE statistics, labeled 
“BE” (black lines). In the latter case, the particles have the 
mass meff = ^7 g. 

Whenever the system is in local thermal equilibrium, 
the temperature dependence may be eliminated to con¬ 
struct the equation of state (EOS) of the GZ plasma 

^ GT . — £gz {Pgz) ■ ( 11 ) 

This relation is further discussed in Sec. IIIB. In this 
case, the entropy density can be also calculated from the 
thermodynamic identity 


rji 5 V / 

^ GZ 

which is consistent with the conventional thermodynamic 
identities 

dsoz = Eb^dscz , dJ^Gz = Sczd^G^ . (13) 

In the opposite case, when the system is not in local 
equilibrium, one introduces an effective temperature of 
the system by the Landau matching condition 

£Gz(T)=e. (14) 

In this case, the quantity T is nothing else but a measure 
of the local, non-equilibrium energy density. 

The thermodynamics of the system is governed by the 
IR scale 7 g. In the temperature range studied at RHIC 
and the LHC a systematic solution of the gap equation 
suggests 7 g — const [58], so that the energy density is 
not simply proportional to [62, 79]. This is crucial for 
reproducing a non-trivial peaked behavior of the energy 
density in this region [58, 59]. As a matter of fact, in [67] 
we used 7 g = 0.7 GeV to describe the lattice data [68] 
in the vicinity of the phase transition, in agreement with 
[58]. We stress that the SU(3) Yang-Mills theory with 
Tc = 260 MeV can be described only qualitatively in the 


temperature range T ^ (1 — 5)rc [58]. This is, neverthe¬ 
less, expected since our simple approach does not include 
any sophisticated resummation scheme in the calculation 
of the gluon properties. For a more thorough study of 
the transition regime, T ~ Tc, one should also keep in 
mind the alterations of the Yang-Mills vacuum structure 
[64, 65]. Finally, for temperatures much higher than the 
ones indicated above, one should keep in mind the rela¬ 
tion of the Gribov scale to the YM magnetic scale that 
arizes from the gap equation [59, 60]. We leave possible 
further improvements of our model for future investiga¬ 
tions. We only note here that the precise value of 7 g 
is not important for our present discussion. After fixing 
7g = 0.7 GeV, the characteristic temperature for phase 
transition in our model is about 100 MeV. 

In Fig. 1 we show the temperature dependence of the 
energy density £gz (red solid line) and pressure Pgz (red 
dashed line) of the GZ plasma. For the temperatures ex¬ 
ceeding the value of 7 g, one observes that the GZ plasma 
behaves similarly to the ideal gas as the energy density 
equals three times the pressure [58]. For comparison, in 
Fig. 1 we have also shown the results for the energy den¬ 
sity (black solid line) and pressure (black dashed line) of 
the massive BE plasma. It is natural to use the effective 
mass 

TOeff = V^7G, (15) 

in the comparative calculations involving the standard 
massive case because it corresponds to the minimum 
of the Gribov dispersion relation for jfej = 7 g, where 
E = •\/2 7g- By inspection of the results presented in 
Fig. 1 we conclude, however, that the use of the effec¬ 
tive mass, given by Eq. (15), is not sufficient to rec¬ 
oncile the results obtained for the GZ plasma and the 
massive BE plasma for moderate and low temperatures. 
The agreement at high temperatures follows from the 
fact that all particles may be treated as massless in this 
case. This is the Stefan-Boltzmann limit where we have 
£ = 3P = 3csbP"^ and the entropy is Ssb = with 

the Stefan-Boltzmann constant Csb = /'^O. Further 

differences are discussed in Secs. IIIB and VG. 


C. Implementation of boost-invariance 

Below we consider a (O-l-l)-dimensional [(O-l-l)D], 
boost-invariant and transversally homogeneous system 
corresponding to the Bjorken model [69]. It is described 
by the flow vector = (t/r, 0, 0, z/r). In this case, all 
scalar functions of space and time depend only on the lon¬ 
gitudinal proper time t, defined as r = \/t'^ — z^. One 
may furthermore introduce the boost-invariant variables 
[80] 

V = kPt — k\\z , (16) 

w = k\\t — k°z . (17) 

The variables v and w replace and /cy — for a boost- 
invariant system one can consider all the quantities in 
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the plane z = 0 where t = t, v = k^t, and w = feyt. We in terms of their respective shear and bulk viscosities 
find that using relations valid at first order, 


k ■u = ^ ^ , (18) 

so that the particle energy becomes a function of the 
boost-invariant variables t, w, and k±, 


E{k-u) = E{T,w,k±) = +fci + ^2^2 1,2 

y ^ ' ^_L 

The phase space integration measure can be written as 

/ ^ /oo ^ 

We note that the phase space distribution function, which 
is a Lorentz scalar, may depend only on r, w and k±, 
namely / = /(t, w, A;_l) [80]. 

For completeness, we also write down the expressions 
for the longitudinal and transverse pressure components 
that will be used in the following [81] 


^11 = 
P± = 


dK 


dK 


t'^E{t, w, k±) 


kl 


2 E{t, w, k±) 


1 - 


1 - 


7^ 


(u;Vr2 + fc^)2_ 


7^ 


(wyT^ + ki)^ 


f, (21) 
f, (22) 


such that the total pressure, P = (Py -I- 2P±)/3, corre¬ 
sponds to Eq. (9). 


III. FIRST AND ZEROTH ORDER 
HYDRODYNAMICS FOR (0+l)D SYSTEMS 

A. Energy-momentum conservation 

The main dynamic equation governing the GZ plasma 
follows from the energy and momentum conservation. 
Using the energy density, e, that is always kept equal 
to its equilibrium value by force of the Landau matching 
condition in Eq. 


see e.g. [82, 83]. 
as 

de 

dr r t2 


(14), for (O-l-l)D systems we can write 
de e -I- P| 


d7 + 


= 0 , 


(23) 


Equation (23) can be further rewritten 


e -|- Pqz 3?7eff + Ceff _ 



VeS 

T 

Ceff 

r 


(27) 

(28) 


where P = (Py -I- 2 P_l)/ 3. For conformal systems the 
shear tensor is traceless, 11 = 0, which immediately 

signals a vanishing bulk viscosity. 

Equation (23) can be treated in two different ways. In 
the first case it is treated as a consequence of the rela¬ 
tivistic Boltzmann kinetic equation which yields e, Py, 
and P_L as functions of the proper time. Once e, Py, and 
P_L, are known one can calculate the equilibrium pressure 
Pgz (if file EOS is known, which is usually assumed) and 
determine the effective coefficients rjes and Ceff as func¬ 
tions of time. 

In the second case, Eq. (23) is considered as one of the 
hydrodynamic equations, given by Eqs. (23)-(28), which 
can be solved only if the kinetic coefficients t] and C are 
known. Such equations are called the first-order (Navier- 
Stokes) fluid dynamic equations. In general, they are 
known to suffer from several deficiencies, but they usu¬ 
ally represent a good description of systems that are very 
close to equilibrium. For such systems one determines the 
kinetic coefficients rj and ( as functions of T. We skip the 
labels “effective” in this second case to stress that a differ¬ 
ent method of determination of the kinetic coefficients is 
used in this situation. Namely, instead of solving the ki¬ 
netic equation for genuine non-equilibrium systems, one 
typically considers small perturbations around the equi¬ 
librium state to determine the linear response of the sys¬ 
tem. This method will be used to determine r] and for 
the GZ plasma in Sec. V. 

The key point of our approach is that the two ways of 
treating the shear and bulk viscosities explained above 
should become equivalent for systems approaching local 
equilibrium. In Section IV we construct the kinetic equa¬ 
tion and find its solutions. This allows us to determine 
rjeS and ^eff as functions of time. On the other hand, 
analyzing small perturbations of the equilibrium distri¬ 
butions in Sec. V we find r] and C as functions of temper¬ 
ature. The consistency check that ?7eff(T) = and 

Ceff('r) = C{T{t)) demonstrates an overall consistency of 
our framework and supports our derivation of the func¬ 
tional forms for C{T) and ri{T). 


B. Bjorken hydrodynamics 


Here we have expressed the shear and bulk viscous pres¬ 
sures [84], 


{P^-P\\) - 

(25) 

~ Pgz 1 

(26) 


For perfect fluid (zeroth-order) hydrodynamics the 
transport coefficients vanish, ( = t] = 0, and the (O-l-l)D 
system is described by the well known Bjorken evolution 
equation [69], 

d£oz(T(r)) £Gz(r(r)) + Poz(r(r)) 
dr 


T 


0. (29) 
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FIG. 2. (Color online) Temperature dependence of the sound 
velocity squared in the GZ plasma (red, solid line) and in the 
massive BE plasma (black, dashed line) with the particle mass 
given by Eq. (15). The green, dotted line represents the low- 
temperature limit Cs =T/meff = T/(%/2 7g), while the red, 
dotted line and the black, dashed-dotted line depict the first- 
order deviation at high-T for the two systems, respectively, 
see Eqs. (35) and (36). 


The solution of this equation will be denoted later as 
Tgz{t). Using the thermodynamic relations in Eqs. (12) 
and (13), we find from Eq. (29) that 


dsr 


dr 


^=0, 

r 


which has the scaling solution 


5gz('^) — 'SGz('n3) 

T 


(30) 


(31) 


This describes the characteristic drop of entropy density 
which is independent of the EOS of the system. The 
time dependence of the temperature in the Bjorken model 
Tgz(t) is found from the expression 


where 


dlnTGz(r) _ 2 

d In r ’ 

(32) 

2 _ dP^z _ 

dscz deczIdTaz 

(33) 


is the speed of sound of the GZ plasma. For confor¬ 
mal systems, e.g., a relativistic, massless gas, where 
e = 3P, we reproduce the well-known scaling solution 
T(r) = To(ro/r)'^= with = i/s- For more general sys¬ 
tems, Eq. (32) demonstrates that the evolution in local 
thermodynamic equilibrium is determined completely by 
the form of Cg. 

The temperature dependence of the sound velocity 
squared for the GZ plasma is shown in Fig. 2 (solid red 
line). The result for the GZ plasma is compared with the 
result for the massive BE plasma (black dashed line). In¬ 
terestingly, the two results differ for moderate and high 


values of T. At low temperatures the common limit for 
the two systems is (see [85] and Appendix G 2) 


cliT) = 


T 

rUeS 


T 


(T ^ 0). 


(34) 


The behavior at high temperatures differs notably. For 
the GZ plasma at T 7 g we have [58] 


2 _ 1 
“ 3 


1 - 




3'\/2^ CgB 
while for the BE plasma we have 


+ 0 


^ ) - (35) 


c 


2 

s 


1 

3 
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(36) 


see Appendix E. The high-T expansion breaks down for 
F ^ 7 g or T < TTieff in the two cases, as seen in Fig. 2. 
The GZ plasma approaches the Stefan-Boltzmann limit 
Cj = 1/3 much faster than the BE plasma, since the devia¬ 
tions go like the third rather than the second power of the 
inverse temperature. The difference between the two sys¬ 
tems under consideration is however mostly pronounced 
in the intermediate temperature regime, specifically for 
F ~ 7 g) where the GZ plasma has a more pronounced 
change of behavior. These features reflect a harder EOS 
of the latter compared to the BE plasma and result in 
a more rapid change of the energy density and pressure 
around the transition temperature, in line with expecta¬ 
tions from the lattice [68]. The full, numerical result for 
the speed of sound is however crucial in order to describe 
the temperature regime relevant for present-day collid¬ 
ers, where both the low- and high-T expansions fail, see 
Fig. 2. 


IV. KINETIC EQUATION IN THE 
RELAXATION-TIME APPROXIMATION 

A. Relaxation time approximation 

Our main starting point for kinetic considerations is 
the formula for the energy density (4), which in the 
(O-l-l)D case depends only on proper time, e(r) = 
/dAT(r, w,k±_)f{T,w,k_\_). In order to study the time 
evolution, we take the derivative of the energy density 
with respect to the proper time, which gives 

de g(T) + T||(t) 
dr r 

= / dK E[T,w,k±) - --■. (37) 

Here we have identified the term proportional to the lon¬ 
gitudinal pressure. We point out that additional terms 
would appear in Eq. (37) if the Gribov parameter 7 g is 
allowed to be medium dependent. 

The terms on the left-hand side of Eq. (37) should 
vanish due to energy-momentum conservation (23) in a 
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(0+l)D system. Thus, the term on the right-hand side of 
Eq. (37) should vanish as well. This suggests that we can 
use the standard kinetic equation in the relaxation-time 
approximation (RTA) of the form [86-88] 


where the damping function D(t 2 ,ti) has the form 


D{t 2 ,ti) = exp 



dr 

Ael('r). 


(41) 


df{T, w, k±) _ /gz(t, w, k±) - f{T, w, k±) 
dr Tredr) 

where 

J dK E{T,w,kj_) fczir, k±) 

= J dK E{T,w,k±) f{T,w,kj_) . (39) 


In Eq. (39) we recognize the Landau matching condition 
for the energy density, see Eq. (14). 

At this point it is important to emphasize that al¬ 
though Eqs. (38) and (39) are sufficient to construct a 
consistent (albeit simple) kinetic theory scheme, they 
are limited to the case of a time-independent relax¬ 
ation time. In particular, Eq. (39) does not hold for 
a momentum-dependent relaxation time. For that to be 
the case, the relaxation time should be included in the 
integrands in Eq. (39). The incorporation of the mo¬ 
mentum dependent relaxation time would make our re¬ 
sults more realistic by implementing the fact that parti¬ 
cles with larger (smaller) momenta interact more weakly 
(strongly). However, such a generalization would make 
our considerations of exact solutions of the kinetic equa¬ 
tion much more complicated (not to speak about the 
fact that at the moment very little is known about the 
scattering of Gribov gluons). Nevertheless, although the 
present approach is restricted mainly to the case of a 
momentum-independent relaxation time, it can be ap¬ 
plied easily in the situations where the relaxation time 
depends on the effective temperature of the system fol¬ 
lowing Refs. [89, 90]. Possible effects of inclusion of the 
momentum dependent relaxation time in our framework 
are discussed in Sec. VD and in Appendices D and E3. 

Another simplification done implicitly in Eqs. (38) and 
(39) is that the relaxation time accounts mainly for in¬ 
elastic collisions. This follows from the assumed form 
of the equilibrium distribution function (6), where the 
chemical potential has been set equal to zero. If the 
system was dominated by elastic collisions, the number 
of particles would be (approximately) conserved and the 
background equilibrium distribution (6) would include a 
non-vanishing chemical potential determined from an ad¬ 
ditional Landau matching condition (similar to Eq. (39) 
but expressing the equality of particle densities rather 
than particle energies). 

The formal solution of Eq. (38) with a time dependent 
relaxation time is [89-91] 

f{T,w,k±) = fo{w,k±)D{T,To) 

+ [ ,, A>('r,r')/Gz(r',w,fc_L), (40) 

Jro Trel[T') 


If the relaxation time is extremely short, the form of 
Eq. (38) guarantees that the actual distribution function 
is in practice always equal to the equilibrium one. 

In order to construct the solution of Eq. (40) we have 
to know the dependence of T and Trei on the proper 
time T. We also have to know the initial distribution 
foiw-,k±), which for this work will be chosen simply 
as the isotropic GZ equilibrium distribution as given in 
Eq. (6). The proper-time dependence of temperature is 
found self-consistently from the Landau matching condi¬ 
tion, see Eq. (43). 

Although the relaxation time Trei in (40) may depend 
on time, in this work we fix it to be a constant in order 
to single out basic features of the modified IR dynamics. 
In the following, we present results for a range of relax¬ 
ation times, 0.5 fm/c < Trei < 2 fm/c, which are in 
the ballpark of phenomenologically relevant values [89]. 
In order to avoid being biased by this choice we have 
also calculated the ratio of the bulk to the shear viscosi¬ 
ties, C/Vj which is independent of the relaxation time and 
contains information about the dynamics of the system 
that is independent of the RTA scheme, see Sec. V. One 
should note, however, that the RTA itself allows for such 
independence. 

The non-equilibrium evolution equations for the energy 
density, the longitudinal, transverse and total pressures, 
I = {e,P||,P_L,P}, respectively, are given by the univer¬ 
sal formula 


I(t) = D{t,to)Hx 
dr' 


7g Tq 
Tiro)' T 


-b 


Trel(T' 





(42) 


where the auxiliary functions ffx are listed in Ap¬ 
pendix B. We choose the initial condition to be Tg = 
'r(To) = 0.6 GeV and tq = 0.5 fm/c, as in [67]. 

As stated above, the temperature of the system is fixed 
by the Landau matching condition which makes sure 
that the actual distribution function f(T,w,k±) yields 
the same energy density as the GZ equilibrium function 
fGz{T,w,kx)- This gives the implicit equation [89-91] 


He 



D{t,to)H^ 



dr' 

Trel(r') 




( 

\T{r,y r) 



(43) 


which has to be solved for the temperature T at the 
proper time r, T{t). Alternatively, in equilibrium, the 
proper-time dependence of the temperature can be found 
directly from Eq. (32), which corresponds to taking 
rrei ^ 0 in Eq. (43). 
















(a) 



FIG. 3. (Color online) The proper-time dependence of the 
effective temperature T (left panel) and the same normalized 
by the equilibrium temperature Tgz (right panel) for three 
values of the relaxation time, Trei = 0.5 fm/c (red, solid), 
1 fm/c (green, dashed) and 2 fm/c (blue, dotted), respectively. 
The equilibrium temperature in the GZ plasma is plotted as 
the dashed-double-dotted line in the left plot, where we also 
have added the temperature dependence for an relativistic, 
massless gas (black, thin, solid line). We will keep these plot¬ 
ting conventions for all subsequent figures. 


B. Numerical results 

In Fig. 3a we plot the proper time dependence of 
the temperature for a massless, relativistic gas (black, 
solid, thin line), for the GZ plasma in equilibrium (black, 
dashed-double-dotted line) and for the GZ plasma with 
three chosen values of the relaxation time, Tj-ei = 0.5, 1 
and 2 fm/c. In the lower panel. Fig. 3b, we plot the ratio 
of the effective temperature to the equilibrium tempera¬ 
ture, T/Tgz, for the same choice of the relaxation times. 
With increasing relaxation time we find slower temper¬ 



T [GeV] 


FIG. 4. (Color online) The interaction measure, Eq. (10), 
scaled by the fourth power of temperature, I/T*, as a function 
of the effective temperature. The conventions for the labelling 
of the curves adopted from Fig. 3. 


ature evolutions and larger T/Tgz ratios. This effect is 
due to increased entropy production within the system 
caused by dissipation. 

In Fig. 4 the interaction measure (10) scaled by the 
fourth power of temperature, 1/T^, is plotted versus T. 
The three curves corresponding to the three different re¬ 
laxation times agree at high and low temperatures. The 
agreement at high temperature follows from using the 
same initial condition in the kinetic equation. The agree¬ 
ment at low temperatures is a consequence of the fact 
that the analyzed systems are all evolving towards local 
equilibrium for sufficiently large proper time. Indeed, to 
cool the system from the initial temperature Tg = 600 
MeV down to T = 200 MeV takes more than 20 fm/c, 
which is much larger than the largest relaxation time used 
in the calculations. Interestingly, the non-equilibrium 
evolution makes the interaction measure always larger 
than the equilibrium case [67]. 

Dynamic properties of the GZ plasma are further de¬ 
scribed in Fig. 5 where we plot the proper time evolution 
of the total, longitudinal and transverse pressures. The 
results are presented again for three values of the relax¬ 
ation time. In all the cases we observe an initial decrease 
of the longitudinal pressure and increase of the transverse 
pressure. Usually, this effect is explained by the presence 
of the shear viscosity that, due to strong longitudinal 
expansion, decreases the longitudinal pressure and in¬ 
creases the transverse pressure [82]. The differences in 
pressure become more pronounced with increasing vis¬ 
cosity, hence with the relaxation time. This is also seen 
in our calculations. We note, however, that the results 
in Fig. 5 represent exact solutions, cf. Eq. (42), of the 
kinetic equation and do not necessarily rely on the inter¬ 
pretation related to the use of first-order dissipative fluid 
dynamics. We also note that the pressure components 
become approximately equal at a relatively late stage of 
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FIG. 5. (Color online) Evolution of the various pressure 
components in time: the total (solid lines), longitudinal (dot¬ 
ted lines) and transverse pressure (dashed lines), respectively. 
The results are presented for three values of the relaxation 
time. The color coding is the same as in Fig. 3. 


the evolution of the system. In contrast, the agreement 
between the full kinetic theory and first-order dissipative 
fluid dynamics sets in relatively fast, approximately at 
T ^ (3 — 5) Trei, See Sec. V. 


V. TRANSPORT COEFFICIENTS 

The basic non-equilibrium properties of the GZ plasma 
at late times are characterized by transport coefficients. 
In the absence of any conserved charges, the most impor¬ 
tant ones are the bulk and shear viscosities, see Eq. (24). 
In this Section we recall our results on the bulk viscosity 
C, obtained in [67], and present new results on the shear 
viscosity t] and the ratio C/v- In nil thn cases we first in¬ 
troduce the effective coefficients using the solutions of the 
kinetic equation (38). Then, small perturbations around 
the equilibrium are studied in order to derive formulas 
for the transport coefficients that define C and rj as func¬ 
tions of the temperature T and the Gribov parameter 
7 g. Numerical results, based on Eq. (43), demonstrate 
the consistency of the two approaches whenever the sys¬ 
tem is close to equilibrium. 


A. Bulk viscosity 

Although the energy density is always equal to the 
equilibrium energy density (due to the Landau match¬ 
ing condition), this does not hold for the pressure. The 
difference between the actuaf and the equilibrium pres¬ 


sures describes the bulk viscous pressure, see Eq. (26). 
Using the definition of the trace anomaly, Eq. (10), and 
the Landau matching condition, Eq. (14), one finds that 

which simply corresponds to (7°^ — 7)/3. Knowing 11 
from the kinetic theory, we may extract the effective bufk 
viscosity using Eq. (28), which gives the effective bufk 
viscosity coefficient, 

Ceff(T) =-rn(r). (45) 

Alternatively, we may seek the solution of the kinetic 
equation (38) in the form 

/ ~ /gz + Sf + ■ ■ ■ , (46) 

where df = —T^eiifcz/dr and the ellipsis represents 
higher order terms. We truncate the series at the lin¬ 
earized level in order to study the late-time behavior of 
the system. In this case. 


sf = - 


E "Trel 


Tt 

dlnrl 
dlnr I 


1 - 


7^ 


/ 

V t 2 


kl) 


/gz (1 + fez) ■ 


(47) 


The second term in the curly brackets can be replaced 
by the speed of sound, given in Eq. (32). Substituting 
f — fez = 5f in Eq. (44) and extracting the bulk viscosity 
from Eq. (45), yields a closed expression, 


. ^ 9o1g Ael 
^ 37r2 T 



3y4 + i_ 


/gz(1 + fez), 


(48) 


where /gz = {exp[jG\/y'^ +y-'^/T] - 1}“^ and 
y = {k ■ u )/ 7 g . We have omitted the subscript “eff” 
since Eq. (48) results from the standard definition of the 
bulk viscosity that applies to first-order hydrodynamics. 

In Fig. 6 we plot the bulk viscosity obtained numer¬ 
ically from Eq. (45) and scaled by the entropy density 
and the relaxation time Trei, C/(s'7-ei), versus T (filled 
symbols). Since we are interested in the properties of the 
system close to equilibrium, we simply approximate the a 
priori non-equilibrium entropy s by the relation [e+P) jT 
(which is strictly only valid in local, thermal equilibrium). 
The results are compared with the linear-response result 
(red solid line), from Eq. (48). We observe a good agree¬ 
ment between the two calculations at low temperature, 
that correspond to large evolution times (much larger 
than the relaxation time). By comparing the results ob¬ 
tained with three different relaxation times, see Fig. 3a, 
we estimate the time when the system complies with a 
description in terms of linearized transport coefficients as 
~ 5 Trei. For all three relaxation times, the ratio C/(sTrei) 
is non-monotonic and exhibits a maximum at a similar 
temperature as the maximum of the interaction measure, 
cf. Fig 4. 
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FIG. 6. (Color online) Effective bulk viscosity scaled by 
entropy density and the relaxation time as a fnnction of the 
temperature for three relaxation times: Trei = 0.5 fm/c (red 
squares), Trei = 1 fm/c (green circles) and Trei =2 fm/c (blue 
triangles). Analytic results in the linearized approximations 
(red solid line) are calculated using Eq. (48). 


It is also instructive to determine the high-T depen¬ 
dence of the ratio C/s. In order to find the leading be¬ 
havior it is simply sufficient to substitute the entropy 
density and the speed of sound by their respective values 
in the Stefan-Boltzmann limit, i.e. s « Ssb and « 1 / 3 . 
In this case, i.e., for T ^ 7 g (with constant 7 g and Trei), 
we obtain 


c 

s 


5 7G'^rel 
8V2tt3 


+ 0 



(49) 


B. Shear viscosity 

The shear viscous pressure arises in locally anisotropic 
systems. In (O-l-l)D systems, this is measured by tt, de¬ 
fined by Eq. (25). From Eq. (27), we find that 


FIG. 7. (Color online) The same as Fig. 6 but for the shear 
viscosity. 


where the definitions for /gz and y are given below 
Eq. (48). Equation (51) is one of the main new results 
presented in this work. 

In Fig. 7 we plot the shear viscosity to entropy density 
ratio, where 77 was obtained numerically from Eq. (50) 
and compare it with the result obtained using Eq. (51). 
We employ the same approximation for the entropy den¬ 
sity as in the previous subsection. The presentation of 
the results adopts the labeling conventions of Fig. 6 . We 
observe a very good agreement between the two calcu¬ 
lations at proper times much larger than the relaxation 
time. While the linearized regime sets in slightly faster 
for the shear viscosity coefficient than for the bulk vis¬ 
cosity one, roughly at 3Trei, there is an overall similar 
systematics for the two quantities. 

In order to find the high-T limit for the shear viscosity, 
it is useful to rewrite Eq. (51) as 


-q = 


50TrelT4„ /7 g\ 

307r2 \t) ^ 


(52) 


= =-'^{P\\-Pa.) ■ (50) 

Whenever the system is close to equilibrium, we may 
again analyze small perturbations around the equilibrium 
distribution, which are described by Eqs. (46) and (47). 
The leading-order term /gz is isotropic and does not con¬ 
tribute to the difference Py — Pj_. Similarly, the term 
d In T/d In T in 6f is also isotropic and can be neglected. 
The remaining integral defining q can be written in the 
form 


where 


/C(a) 



+ a* (e^ - 1)^ ’ 


(53) 


where v = ^gU/T and A = -\- a'^ jv^. For asymp¬ 

totic temperatures we simply pick up limT->.oo f] = 
/C(0), where /C(0) = 47r"‘/15. Thus, since in the 
same limit the entropy of the GZ plasma simply goes to 
the Stefan-Boltzmann limit, we get that 


q = 


9 o 1 g 'Tel 
307r2 T 


' (/-I)' 

/gz(1 + /gz), (51) 

2/^ + 1 


7 _ TrelT 

S “ 5 


+ 0 



(54) 
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FIG. 8. (Color online) The ratio of the bulk and shear 
viscosities, scaled by 1/3 — for the GZ plasma (red, 
solid line) and the BE plasma (black, dashed line). The black, 
horizontal line indicates the value 5/2. 


FIG. 9. (Color online) The temperature dependence of the 
bulk to shear viscosity ratio n x (^/rj for the GZ plasma, see 
Eqs. (48) and (51) with n = 1 (red, solid line), and the massive 
BE plasma, see Eqs. (E7) and (E8) with n = | (black, dashed 
line). 


For a massless, relativistic ideal gas, this relation is exact 
at any temperature in the RTA [89]. Naturally, this limit 
also coincides with the result for the massive BE plasma, 
see Sec. E. 


C. scaling 

Erom Eqs. (48) and (51), which were obtained assum¬ 
ing a constant Gribov parameter and relaxation time, we 
conclude that the ratio (jt] is independent of the choice 
of the relaxation time. As a consequence, this ratio is 
a general prediction of our approach. Strikingly, while 
typical perturbative arguments predict oc (Cg — 1 / 3 )^ 
[12, 92, 93], we obtain a different form of the scaling. 
This is most easily demonstrated in the high tempera¬ 
ture limit, r 7 g, where, combining Eqs. (49) and (54) 
which result in (//?7 « (^) , with Eq. (35), gives 

^ = Kgz Q - -b ... (T » 7 g) , (55) 

with Kgz = 5/2 and the ellipses stand for terms that are 
power-suppressed in the ratio 7 g/ 7’- In Fig. 8 we plot C/?? 
for the GZ plasma using Eqs. (48) and (51) (see the solid, 
red curve), where the limiting behavior is clearly seen at 
high temperature. The linear scaling (55) is characteris¬ 
tic for systems where the conformal symmetry is explic¬ 
itly broken and has also been found in strongly coupled 
theories [44, 45]. As an interesting observation, we note 
that the agreement between the Gribov-Zwanziger quan¬ 
tization at high temperature and strongly-coupled theo¬ 
ries based on gauge-gravity duality is not a singular case, 
as the recently discovered massless mode of the QGP [61] 
is also closely in line with the holographic quasinormal 
mode [94]. 


It is instructive to contrast Eq. (55) with results for 
the BE plasma. The bulk viscosity of the BE plasma 
scales very similarly with I — Cg as the result for the GZ 
plasma, see Eq. (E7) and the shear viscosity behaves as 
in the ideal gas at high temperature, see Eq. (49). This 
is in agreement with the results found in [95], which has 
also been derived within a quasiparticle model using RTA 
[96]. However, we find that the high-T limit, T Weff, 
results in ^ « yUa while the speed of sound 

deviates from its ideal value by a quadratic term, see 
Eq. (36). Therefore, we conclude that for the BE plasma 
at high temperatures, we find the relation 

c _ 

— ^BB 

v 

with Kbe = 31 / 15/2 ~ 5.81, which is qualitatively differ¬ 
ent from the linear relation found for the GZ plasma. The 
difference to the expected scaling with the second order 
[12, 92, 93] can be traced to the lack of a momentum- 
dependent relaxation time as well as a temperature- 
dependent quasi-particle mass, as expected from pertur¬ 
bative arguments. We plot (//yy for the BE plasma using 
Eqs. (E7) and (E8) (see the dashed, black curve) in Fig. 8. 
We note that the curve monotonically tends to zero at 
high temperature, confirming the absence of linear scal¬ 
ing. 

To further study the differences between the behavior 
of the GZ plasma and the massive BE plasma, in Eig. 9 
we plot the //ry ratio versus T for these two systems. 
At low temperatures, T ^ 7 g or T ^ to, the C,/ri ratio 
tends to 5/3 for the GZ plasma and to 2/3 for the massive 
plasma, see Appendix G 5 for the results of an analytic 
low-T expansion. In terms of the constants in the respec¬ 
tive scaling behaviors, this corresponds to Kqz = 5 and 


3 


0/ Z 


+ . . . (T > TOeff) , (56) 
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Kbe = 2, since —)• 0 in both cases. For a comparison of 
their respective low-T behaviors, we have multiplied the 
BE result by 5/2 in Fig. 9. Surprisingly, the upper bound 
of C/j? in the GZ plasma is in good agreement with the 
maximal value of this ratio from a recent gauge-gravity 
study [47]. 

In order to access the phenomenologically relevant 
temperature regime we have to use the full expressions, 
Eqs. (48) and (51). In this regime, T ~ (1 — 5)T'c, the 
scaling relation in Eq. (55) is only approximate with a 
slowly varying see Fig. 8 . We note, however, that 
the C/r; ratio in the GZ plasma is always significantly 
larger than in the BE plasma. In particular, at small 
temperature C > r] in the GZ plasma while this is never 
the case in the massive BE model. 


D. Comments on the momentum dependence of 
the relaxation time 


Our main results have been obtained with two impor¬ 
tant assumptions, namely, that the Gribov parameter 
and the relaxation time are treated as constants. One 
may ask the question how these two assumptions affect 
our considerations. In order to analyze this issue, in Ap¬ 
pendix D we relax the latter requirement presenting cal¬ 
culations done with a momentum-dependent relaxation 
time. In this case, we assume a simple power-law rela¬ 
tion 


T-rel = T-°el(|fc|/7G)“, (57) 

with a > 0 being a parameter (for other discussions of the 
momentum dependent relaxation time see, for example. 
Refs. [97, 98]). 

The relaxation time (57) should appear in the Landau 
matching condition (39) that had been used before to ob¬ 
tain Eqs. (44) and (48). With a more complex form of 
the relaxation time, we have to use a different strategy 
and employ (47) directly to calculate first the bulk vis¬ 
cous pressure and then to obtain the bulk viscosity. In 
this way, instead of Eq. (48), one finds the expression 


c = 

X 


g07G ^r°el 

67r2 T 
'1 -1 
.32/4 + 1 



d2/2/“ (/ 


1 ) 



/gz(1 + /gz)- 


(58) 


For the shear viscosity we can still use Eq. (51), pro¬ 
vided the relaxation time Trei is moved into the integrand, 
which gives 


go7G A°ei 
' 307r2 T 



d 2 / 2 /“ 


2/4 + 1 


/gz(1 + /gz). (59) 


We note that using the definition of the sound velocity 
one can check that in the case a = 0 Eq. (58) is equivalent 
to (48) with Trel = r°ei. 

Using (58) and (59) we construct the C/g ratio and 
find its asymptotic behaviour in the limit 7 g T. The 


analytic calculations presented in Appendix D show that 
for a > 1 one gets 



For 0 < a < 1 one gets also a linear scaling with the 
factor 5 decreasing to 5/2 as a decreases from 1 to 0, 
which agrees with our previous result (55). Thus we find 
that the scaling of the C/r/ ratio with the first power of 
the conformal measure turns out to be a universal feature 
of the GZ plasma. 

In a recent paper, Jaiswal, Friman and Redlich showed 
that linear scaling holds also for a classical massive gas if 
the relaxation time depends linearly on the particle en¬ 
ergy [98]. In Appendix E we argue that our approach 
leads to identical results if applied to the same system. 
We find linear scaling also for the massive Bose-Einstein 
gas if the relaxation time depends on the momentum ac¬ 
cording to Eq. (E12) with a > 0. However, we note that 
these results have been obtained assuming temperature 
independent Gribov and mass parameters, in the case of 
the GZ and BE systems, respectively, which may affect 
the final form of the scaling. We plan to analyze this 
issue in a separate work. 


VI. SUMMARY AND CONCLUSIONS 

In this work we have continued our earlier investiga¬ 
tions of the equilibrium and non-equilibrium properties of 
a plasma of gluons resulting from the Gribov-Zwanziger 
quantization (GZ plasma) in the relaxation time approx¬ 
imation. We have supplemented our earlier results for 
the bulk viscosity by the calculation of the shear viscos¬ 
ity and the C^jt] ratio as a function of temperature. We 
have demonstrated that the exact solutions of the kinetic 
equation support our expressions for the two viscosities. 

Several of the studied features suggest that at RHIG 
and LHG energies the GZ plasma is a strongly coupled 
system. Firstly, we observe a larger bulk viscosity, 
compared to conventional quasiparticle approaches espe¬ 
cially in the low temperature regime, where C > V- This 
highlights its signihcant role for heavy-ion phenomenol¬ 
ogy. Secondly, for constant 7 ^ and Trei, we find a linear 
scaling of the C/r/ ratio with 1/3 — c^, reflecting strong 
breaking of conformal symmetry of the system. At high 
temperature, this is at variance with the expected behav¬ 
ior of weakly-coupled plasmas [36] and similar to expec¬ 
tations from a strongly-coupled system [44, 45]. 

Our results have been supplemented by the discussion 
of several aspects of the Gribov approach. In particu¬ 
lar, we have analyzed the Lorentz covariance and found 
the low- and high-temperature expansions of thermody¬ 
namic functions and transport coefficients. We have also 
performed systematic comparisons between the qualita¬ 
tively different properties of the Gribov-Zwanziger and 
the massive Bose-Einstein plasmas. In this exploratory 
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study, our results have been obtained under the condi¬ 
tion that the Gribov scale as well as the relaxation time 
are kept constant. Relaxing these assumptions allows to 
study more complex situations. 

Furthermore, we have checked the effects of a 
momentum-dependent relaxation time on our results and 
confirmed a linear scaling of the C /rj ratio with the con¬ 
formal measure also in this case. Similar scalings hold 
also for the massive Bose-Einstein and Boltzmann gas if 
the relaxation time is proportional to the momentum or 
energy of particles. This issue deserves further studies to 
clarify the interpretation of a linear scaling as a signal of 
strongly coupled systems. 

In conclusion, we have demonstrated that the Gribov- 
Zwanziger quantization can be used to address in- and 
out-of-equilibrium physics in a unified way, which can be 
useful for future phenomenological applications to ultra- 
relativistic heavy-ion collisions. The latter may require 
the inclusion of a temperature-dependence of the Gribov 
scale. The first step in this direction has been done in 
Ref. [99]. 
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Appendix A: Implementation of Lorentz covariance 

Gribov’s dispersion relation, Eq. (1), relating the par¬ 
ticle three-momentum k and energy E explicitly breaks 
Lorentz invariance. The reason for this breaking is the 
use of the Coulomb gauge V • A = 0 in the derivation of 
Eq. (1). Without the lack of generality, we may assume 
that the Coulomb gauge is imposed in the inertial refer¬ 
ence frame S. The four-velocity of the whole system in 


S IS = (1,0,0,0), hence, we rewrite Eq. (1) as two 
combined equations 

= u'^ = (l,0,0,0). (Al) 

k 


We may try to rewrite Eq. (Al) as a single covariant for¬ 
mula. The first natural choice seems to be the introduc¬ 
tion of the four-vector k = {Ek,k) with the Minkowski 
square = E^ — k^, then Eq. (Al) can be written co- 
variantly as 


(fc • u)^ — kP' 


(A2) 


In the frame where = (1,0,0,0), Eq. (A2) is reduced 
naturally to (Al). However, when (A2) is considered in 
other frames, it typically yields four different solutions 
(some of them complex) for the energy Ek determined at 
given values of k and v. This leads to a paradox, since 
in S all energy states are real and equal. 

Clearly, one should restrict oneself to the energies de¬ 
fined in the local rest frame of the fluid element. A nat¬ 
ural way to achieve this and implement the covariant 
description is to consider k as the spatial part of the four 
vector k = {ko = \k\,k) and to treat the energy E defined 
in (1) as a scalar 


In the frame S, where = (1,0, 0,0), Eq. (A3) is re¬ 
duced directly to (Al). In the frame S' we have 


where u' and k' are the four-velocity of the frame S and 
the four-momentum k seen in the frame S'. For one 
dimensional systems considered in this work, the com¬ 
ponents of u' may be treated as the parameters of the 
Lorentz transformation k' ■ u' = k'^u'^ — k'^u'^ = ko = jfcj. 
Hence, Eq. (A3) in frames other than S gives real and 
non-degenarate values of the energy, in agreement with 
the original definition of the energy done in the frame S, 
where the Coulomb gauge was imposed. 

The prescription (A3) has been adopted in this work. 
We note that our treatment of energy is similar to the 
treatment of temperature in the special theory of rela¬ 
tivity. Temperature measures average kinetic energy and 
has been treated by many authors as the time compo¬ 
nent of the four-vector. However, since T is measured 
in the LRF, its treatment as a scalar quantity seems to 
be more appropriate. The energies appearing in Gribov’s 
formula have been obtained also in the special reference 
frame where the Coulomb gauge condition holds. This 
suggests that E defined in such a way may be treated as 
a scalar quantity in the similar way as T. 
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Appendix B: Auxiliary functions 


The energy density is calculated via the auxiliary func¬ 
tion iJe(7G, a, 5), defined by 




yheiy,b) 


(Bl) 


with 


heiy,b) = b dyj sin ^ (B2) 

Jo P 

where /3^ = 5^ cos^ ip + sin^ p, such that 

h,(y,l) = x/TTT. (B3) 

The energy density in equilibrium is given simply by 

e'^-(T)=i?,(7G/T,l). (B4) 

Next, we work out the auxiliary functions for the longi¬ 
tudinal and transverse pressures Py and Pj_, respectively. 
They are from the expressions 


4 poo 


„ . , ^ 1 doll 


4 poo 


dy 


y fep|| (y, 

^y/y^+V~^ _ 1 ’ 

y (y, 

ayjy^+y-^ _ I ’ 


(B5) 

(B6) 


where 


o sin(/JCOS^ (y^/3^ — l) 

hp„ (y, h) = b^ / dtp- - - , ^ , B7 

’ Jo /33 +1 ^ ^ 

, , , sinV(y4^4_i) 

^P.L (y, b) = b dp ^^ , 

Jo /33vW^+l 


such that 


u ( 1 y"* -1 

«Pii (y. 1) = o 


hp± (y> 1) = 


3 y/y^ + 1 ’ 

2 y"-l 

3 \//n 


(B8) 

(B9) 

(BIO) 


The equilibrium pressure components are then 
P||'=^(T) = Pp,,( 7 g/T, 1) and P_L=^(r) = Ppjyo/T,!), 
and it is straightforward to check that in this case 
P|I<=^(T) = 

The total pressure is given by P = |(P|| -|-2Pj_), where 
the auxiliary function responsible for its evolution is eas¬ 
ily deduced to be 


Hp{a,b) = 


9olt [°° 

67r2 ./n 


dy- 


yhp{y,b) 


(Bll) 


with 


, , sin(^(y4/34_i) 

hp{y,b) = b dp 7^; / , B12 

Jo P^/yW^ 


such that 


hpiyA) = / , ^ ■ (B13) 

\/y^ + i 

The total pressure in equilibrium is then 

P^-{T)=Hp{j^/T,l). (B14) 


Appendix C: Low temperature expansions 

In this Section we analyze different physical quantities, 
such as the equilibrium energy density, equilibrium pres¬ 
sure, the shear and bulk viscosity coefficients, and derive 
analytic expressions for all of them in a form of the se¬ 
ries expansion with the leading terms corresponding to 
the low-temperature limit a = 7g/P 1- By including 
in the systematic way more terms in the expansions, we 
may obtain a successful description also in the region of 
moderate temperatures. 


1. Equilibrium energy density 

We start with the equilibrium energy density. Combin¬ 
ing our earlier results (B1)-(B4) for the case b = (3 = I, 
we find 


^ go7G y\/y'^ + 

27r2 Jo y ^ay/y^+y-^ _ ^ ' 


(Cl) 


By splitting this integral into two parts, the first from 
0 to 1 and the second from 1 to oo, and changing the 
integration variable y to 1/y in the first part we get 


e 


GZ 


go7G dy (y^ -by ^) y/y^ -by"^ 
27 r 2 y ^ay/y^+y-^ _ ^ 


(C2) 


In the next step we introduce the integration variable z 
defined through the relations 


y = \l , 
dy z dz 

y Vz4 - 1 ’ 


(C3) 


Then, we can rewrite Eq. (C2) as 


€ 


GZ 


go7G , fsjz) 

2^2 eV2az _ 1 ’ 


(C4) 


where the function /e(z) is defined through the formula 


fe{z) = {y^ + y ^), (C5) 

V-J — 1 

where y is the function of z defined in (C3). For a 3> 1, 
the main contribution to the integral (C4) comes from 















































15 


the region 2 ^ 1 , where feiz) has the expansion 

, , , 72 237^^ 307(2 - 1)3/2 

f,{z) = + ■ ^ 


272 


1672 


6472 102472 

We represent this expansion in a more general form as 




(C7) 


n =0 


where the coefficients ai" me be read off from (C 6 ). In 
this way we find the expression for the equilibrium energy 
density as the following series 


e = 


9o1g 

27r2 


E 

n—O 


,C 8 ) 


In the next step it is convenient to change the integration 
variable 2 to the variable a defined by equation 2 = 1 + 
c?. In this way we get 


„ -A °° foo In 

.GZ _ 5o7g (n) j ^ 

+ J g+2a(l++) _ 1 

n=0 e ± 


£: = 


TT^ 


(C9) 


This integral over a can be done analytically and ex¬ 
pressed by the gamma and polylogarithm functions, 

_GZ _ 5o7g \-("+5) 


e = 


27r^ 


E4”Mv^« 


n =0 


X r 1 n + - ) Li „+1 




(CIO) 


At very low temperatures, a = Jg/T 3> 1, the most 
important contribution is that with n = 0. In this case 
we get 


go7^ 

(72 a) 


(ClI) 


2. Equilibrium pressure and sound velocity 

Starting from Eqs. (B1I)-(BI4) and following the same 
steps as in the case of the equilibrium energy density one 
finds 


4 poo 4 1 

pGz ^ ^ f g -1 


1 


67r2 


7g^ + y~^ _ 1 


and 


4 /*oo 

pGZ ^ ^ 


GZ g07G / /p(^) 


, (C12) 


(C13) 


67 r2 g+2 az _ I ’ 

where the function fp{z) is defined through the formula 

y5 _^y-5 _ 


/p7) = 


72 72^ - 1 


(C14) 


In the region 2 7 1, fp[z) has the series expansion 
fp{z) = 672 7 ^ + 

93(2-1)3/2 ^ 51(2-1)7/2 


872 


which defines the coefficients Op used in the formula 


3272 

7") 


(C15) 


4 cxo 

pGz ^ go7G ^ 


67r2 


d 2 


_ 1 \n-l/2 


( 2 - 1 ) 


n—1 


c >\/2 az _ 


(CI 6 ) 


Note that Up^ = 0. The integration over 2 can be done 
analytically in the same way as in the case of the en¬ 
ergy density, and we find an analogous expression to 
Eq. (CIO), namely 


4 00 

pGZ = & y 

67r2 21 ^ ^ 

n—0 


X r I n + - ) Li„pi 


2 a 


, -("+ 3 ) 




(C17) 


At very low temperatures, a = ^g/T ^ 1, the leading 
term in (C17) is that with n = 1, 






= -%/ 2 c 


72 7r3/2 


(C18) 


Equations (CIO) and (C17) can be used to determine 
the sound velocity in the system. In the case of the tem¬ 
perature independent Gribov parameter 7 g we obtain 


7 = 


dpGz dP^’^jda 


QgGZ 


de^z j do 


(CI9) 


At low temperatures we may use (Cll) and (C18) as the 
approximations for (CIO) and (C17). Keeping in mind 
that a = 7 g /T, we find 


cliT) = 


(C20) 


1 _ 1 T 

72 a 72 7g ’ 

We note that this result is similar to the result obtained 
for the massive plasma where we find = T/m [85] , 
where m is the mass of particles in the plasma. 

One may go beyond the leading-order term in the 
expansions (CIO) and (C17) and include higher-order 
terms in n in order to describe the regime close to the 
phase transition. Taking Li„_|_i 
to n = 4, we get 

„ — %/ 2 a 




0 -+ 2 C 


up 


3o7g 




1 +? 


8 V72( 


921 


1 

72 a 
175035 / 1 


+ 


11085 / 1 
1024 (vTfa 


32768 \72a 


(C21) 
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FIG. 10. (Color online) The a = 'JgIT dependence of the 
exact speed of sound squared compared to approximated for¬ 
mula (C23) up to n = 1 (green dashed line) and n = 4 (blue 
dotted line). For completeness we show corresponding re¬ 
sults for n = 1 (purple dashed-dotted line) and n = 4 (black 
dashed-double-dotted line) obtained by including first NLO 
term in Li„+i 


and 


pc 


5o7g 


o--\/2 a 


V2 7r3/2 


3/2 


' 19/1 

1 + TT 


8 \y/2i 


465 


1 


128 \^/2a 
and accordingly 


1785 / 1 


1024 Vy2a 


ciiT) = 


1 + 


1 

\f2a 
29/1 


2 ^72^ 

63 


8 \i/ 2 a 


8 Vy 2 c 


(C22) 


(C23) 


In Fig. 10 we present the comparison of the approximated 
formula (C23) with the exact one. We observe that, as 
expected, the expansion (C23) works very well already at 
leading order in n (dashed green line) if the temperatures 
are very small. However, at the same time, we see that 
the expansion completely fails at temperatures close to 
and above the phase transition (a 3), even when in¬ 
cluding higher n terms (dotted blue line). In fact, one 
can show that in order to obtain successful description 
at the phenomenologically relevant temperatures which 
are achievable at RHIC and LHC, 1 ^ a ^ 3 (yellow 
shaded band in Fig. 10), one has to include NLO terms 

in Li„^i expansion. In Fig. 10 we present the 

expansion analogous to (C23) obtained assuming that 
Li„^i -I- e“ 2 v^a 2 -(n-i-i/ 2 ) taking 

terms up to n = 1 (purple dashed-dotted line) and n = 4 


(black dashed-double-dotted line) in (CIO) and (C17). 
We see that see that taking already first NLO term in 
Li„^i is enough for proper description of tem¬ 

perature dependence of the speed of sound close to phase 
transition. We do not show these, somewhat lengthy, ex¬ 
pressions for herein. 


3. Shear viscosity 


For the shear viscosity coefficient we use the for¬ 
mula (51) 


77 = 


1 gOTc'^rel 
10 37r2 T 


dy 




+ l (e^ — 1)2 


, (C24) 


where we have introduced the notation A = a a/ y'^ + y 
This may be further rewritten as 


77 = 


1 901 g '^rel 
10 37r2 T 

with the dimensionless quantity 


77 


’'“X (.^-1) 


e 

T _ 1'|2 


(C25) 


(C26) 


Changing first the 77 -integration range to the interval 
1 < 7/ < 00 and, second, changing the integration vari¬ 
able from y to z using the formulas (C3) we find 


Here 


77 = 


fviz) = 


dzf^iz)- 




^^y/2az _ 

zy (y^ - 1)^(1+ 2/"®) 


— 1 


y^ + 1 


(C27) 


(C28) 


with the expansion 

00 

/,(z) = £a(")(z-l)"-i/^ 

n—1 

= + 34(z - 1)3/2 


111, ^, 5/2 257, ^, 7/2 


(C29) 


In this way we come to the series of the integrals of the 
form 


77 = V 


V 

n—1 


dz 


_ ^'^n— 1 / 2 az 


(z-l) 


-E 


n—1 


(271-1), 
•\/2 a 


^in) .00 ^ 2 („-l) 

- / da 


gV2a(l+a2) _ ^ ■ 


(C30) 


The second line in (C30) was obtained by the substitution 
z = 1 -b a^ in the first line and by integration by parts. 
The integration over a in (C30) yields the gamma and 
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polylogarithm functions, hence, the final result for the 
shear viscosity may be written as the following series 


V = 


1 goTc •^rei 
10 37r2 T 


X r n + 


1 


OO 

n—1 

'■-4 ■ 




(C31) 


To get (C31) we used (2n — l)r(n — 1/2) = 2r(n + 1/2). 


4. Bulk viscosity 


The bulk viscosity coefficient has been given in (48) 
and can be written as the sum of two terms 


> _ SobG/tel / 2 . _ G'l 

^ 37r2 T 3 ' ■ 


(C32) 


The coefficients Q {i = 1, 2) are given by the integrals 



FIG. 11. (Color online) The a = ^gIT dependence of 
“C/h/(Cs —1/3) at leading order (solid red line) and assuming 
Li„+i + e- 2 v^“ 2 -("+i/ 2 ) dashed- 

double-dotted line), see discussion of Fig. 10. 


0 



dz /c. (z) 


oV2c 


j\/2 az _ 


1 )^ 


(C33) 


5. (^/r) ratio for T —> 0 


where the functions and are defined by the fol¬ 
lowing expressions 




(C34) 


n—0 

1 


1 - 1/2 


. + —^- (z - 1)'^/"' 

4 32 ^ ^ 

+ —(z-l)'/2-—(z-1 )^/2+ ■ 
128^ ’ 2048^ 


/c 2 ( 2 :) 


zy ^) 

\/z^ - 1 + 1 


(C35) 


n—0 

! 1 \5/2 / 1 \7/2 


Note that the series defining starts with the term 

2\/z — 1, thus the coefficient vanishes. Collecting 

the expressions for Ci and C 2 together, we arrive at the 
formula 


c 


c; OO 

goTo 'Trel 

37r2 T ^ 



1 

3 




(C36) 


For large a all the polylogarithmic functions can be ap¬ 
proximated by their arguments (this corresponds to the 
limit of Boltzmann statistics). Then, the leading contri¬ 
bution in (C36) is obtained by inclusion of the first two 
terms 


c 


goyc^Treie \ / 1 , _ 1 qA 

37 r 2 T a Vv^a 3 ) 

\V 2 a 4372 y 2 a 


^2 goy^Trei e 

15 ' 


(C37) 


Similar expansion for the shear viscosity (with the term 
n = 1 alone) gives 


_ goTc'^rel e '/2 a 

9 

which leads to the ratio 

C ^ 5 
77 3 ■ 


(C38) 


(C39) 


By comparing this with the scaling relation in Eq. (55) 
we conclude that Kqz = 5. 

Similarly as in Sec. C2 one can go beyond the lead¬ 
ing order in polylogarithm expansion to see how the re¬ 
sult for K presented above is modified when one includes 
higher order therms in a expansion. In Fig. 11 we present 
the K = —C/v/icl — 1/3) at leading order (solid red 
line) of Li„_|_i and assuming Li„_|_i « 
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g-%/ 2 a _|_ g- 2 v^a 2 -(n+i/ 2 ) (bjack dashed-double-dotted 
line). One can observe that close to the phase transition 
temperature the value of k is significantly smaller than 
the LO result and somewhat closer to the value obtained 
in Fig. 8. 


e = +1 for bosons and e = — 1 for fermions. The quan¬ 
tity go is the degeneracy factor connected with internal 
quantum numbers, and g is the chemical potential (set 
equal to zero in the calculations). The series form fol¬ 
lows from the expansion of the BE or Fermi-Dirac (FD) 
distributions, Kn are the modified Bessel functions. 


Appendix D: Momentum dependent relaxation time 


In this Section we give more details describing cal¬ 
culations of the kinetic coefficients of the GZ plasma 
in the case where the relaxation time depends on the 
momentum according to Eq. (57). To find a high- 
temperature leading-order term of the bulk viscosity de¬ 
fined by Eq. (58) we first write 


^ 67r2T 

where 

j/a(y4 

(e^ 


C = / dg 
^0 


3 3(g4 +1) 


(Dl) 


1 ) 

I)^- 

(D2) 
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P 

s 


90 


OO ^ 

90 rj .2 2 Y- e iL^ (m \ 

K—1 

= [(4T - g k) K2 (^k) 

+ itikKi . 


(El) 

(E2) 


(E3) 


It is also worth recalling the high temperature behavior 
of these thermodynamic quantities. At T ^ to we find 
that 


Here we have introduced the notation A = a \/+ y~^ 
and a = 'JgIT. In the high temperature limit, T —>■ oo, 
corresponding to a —>■ 0, we use the approximation 


1 2_ PgzO? _ 15 


(D3) 


Substituting Eq. (D3) into Eq. (D2) one finds that the 
leading term is given by the integral over the variable y 
which is analytic and yields (for a > 1) 


C= ^a-“-2(4 + a)r(4 + a)C(4 + a). (D4) 

Here r(a;) is the Euler F function, while ({x) is the Rie- 
mann () function. In an analogous way, in the case of the 
shear viscosity, we obtain 


rj = 


go7GA°ei - 
SOtt^T 


(D5) 


where the leading term for fj is (again for a > 1) 


f] = a ^ “(4 + a)r(4-I-a)(^(4-I-a). (D6) 


Collecting the results given above we find the scaling 
with the first power of the conformal measure given by 
Eq. (60). We have also found the numerical evidence that 
the linear scaling holds for 0 < a < 1. 


Appendix E: Massive BE plasma 

1. Thermodynamic properties 

In this Section we collect the basic thermodynamic 
formulas for the massive plasma. Below, the parameter 


e = 3csBT4-^m2T2 + 0(l) , (E4) 

P = CsbT 4-^TO^r^+O)!) , (E5) 

s = 4csbT3-^to2t + 0(1) , (E6 ) 

where Csb = 87 r^/45 for SU(3), which reflect the soft EOS 
of the massive BE plasma. The first terms in the above 
expansion are the Stefan-Boltzmann energy density, pres¬ 
sure and entropy, respectively. 


2. Kinetic coefficients 


For vanishing chemical potential g = 0, following the 
same linearization procedure as above, we obtain the bulk 
and shear viscosities as 


> _ gpTO^ Trel 

^ “ 67r2 T 

go TO® Ael 

307r2 T 



1 y'^ 

3g2 + i_ 


/(! + /), (E7) 



g 2 1 


/(! + /), 


(E8) 


where / = [exp(TO-\/g^~+T/r) — 1]“^, y = \k\/m, and 
the speed of sound = 9P/9e|^=o is found from Eqs. 
(El) and (E2), respectively. The viscosities, given by 
Eqs. (E7) and (E8), agree with the ones found from [95]. 

At high temperatures we can also find the generic be¬ 
havior of C/s and g/s. In order to find the leading coeffi¬ 
cient of the high-T expansion, it is enough to simply con¬ 
sider the Stefan-Boltzmann entropy s « SgB = 4 csbT^. 
Then it is straightforward to find the following expres- 
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sions at T ^ m, 


Substituting Eq. (E15) into Eq. (E14) one finds 


c 

s 

1 

s 


5 

167r3 T2 




(E9) 

(ElO) 


which results in the ratio 

C _ 25 /m\^ 
rj V T / 


(Ell) 


at high temperatures. This leads directly to Eq. (56). 


3. Momentum dependent relaxation time 


We have also studied the BE massive gas with the mo¬ 
mentum dependent relaxation time, in the similar way 
as we did it for the GZ plasma in Sec. D. In this case 
{a > 0 ) 

T-rei = r,°i(|fc|/m)“. (E12) 


Following the same steps as in Sec. D we find 


gow^T°i ~ 


(E13) 


where 


c = 


dy 


-c?- 


1 


^ _ ('E14) 

3(2/2 + l)J (e^-1)2-^ 


In this case A = a\J+ 1 and a = m/T. In the high 
temperature limit, T —> oo, we use the approximation 


1 2 _ o _ ^ 


(E15) 


C = [/3BEi^(4 + a) - H{2 + a)], (EI6) 

where H{x) = a;r(x)C(a;). Similarly, in the case of the 
shear viscosity, one obtains 


1 = 


9orrfr^~ 

307r2T 


where the leading term for rj is 


(E17) 


(EI8) 


Consequently, we obtain the scaling 



g(2 + a) W1 _ 2\ 


(EI9) 


As the parameter a grows, the ratio H{2 + a)/H{i + a) 
tends to zero, and the result for the BE massive gas agrees 
with the result for the GZ plasma. On the other hand, 
for a = 0 one finds that H(2)/{P^eH{4)) = 1 and the 
linear scaling (E19) breaks down. In this case one gets 
the scaling of the ratio of the viscosities with the power 
3/2, in agreement with our result given by Eq. (56). 

Finally, we note that we have also done the calculations 
for the massive Boltzmann gas. Interestingly, for a = I 
we reproduce a linear scaling with the prefactor 2, which 
agrees with the recent result obtained in Ref. [98]. 
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